Multi-Species Pair Annihilation Reactions 
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We consider diffusion-limited reactions Ai + Aj — * (1 < i < j < q) in d space dimensions. For 
q > 2 and d > 2 we argue that the asymptotic density decay for such mutual annihilation processes 
with equal rates and initial densities is the same as for single-species pair annihilation A + A —* 0. 
In d = 1, however, particle segregation occurs for all q < oo. The total density decays according to 
a q dependent power law, p(t) ~ t~ a ^ q K Within a simplified version of the model a(q) = (q — l)/2g 
can be determined exactly. Our findings are supported through Monte Carlo simulations. 
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A large variety of systems in physics, chemistry, biol- 
ogy, and ecology can be modeled in terms of diffusion- 
limited reactions. This is because of their unifying fea- 
ture of being composed of mobile agents ('particles') 
which interact upon encounter. The traditional mean- 
field rate equations for such systems apply only to ho- 
mogeneous particle densities. However, in many such 
systems the spatial dimension d has one or more criti- 
cal values below which density fluctuations invalidate the 
rate equations and new phenomena appear. The fluc- 
tuations may stem from, e.g., reaction- induced noise or 
initial state randomness, and typically dominate the sys- 
tem's large-scale, long-time behavior jj]. In order to ex- 
tract the newly emerging nonequilibrium behavior, more 
sophisticated methods are needed: extensive numerical 
simulations (for recent overviews, see Ref. 0) along with 
powerful analytical methods such as scaling approaches; 
mappings to field theories followed by renormalization 
group; and exact solutions (mostly limited to d = 1) ||. 

The theory of annihilation reactions has two land- 
marks, of which one, the single-species pair annihila- 
tion reaction A + A — > 0, represents the simplest model 
case ||. The other one, the two-species annihilation 
A + B — > 0, is considerably more subtle. It exhibits 
the remarkable phenomenon that for d < 4 and for an 
initially random particle distribution with equal densi- 
ties Pa{0) = Pb(0), the two species segregate into pure 
A and pure B domains, and the annihilations become lo- 
calized within sharp reaction fronts between the domains 
H . It is a natural next step to ask for the long-time decay 
properties of a system of q species that mutually annihi- 
late according to Ai + Aj — ► 0, with 1 < i < j < q 
This 'mutual annihilation model' (MAM) is the subject 
of this Letter . A special case is the fully symmetric 
MAM, characterized by equal reaction rates A-y, equal 
diffusion constants Di, and equal initial densities Pi(0). 
We will find that in d > 2 for all 2 < q < oo the fully 
symmetric MAM behaves as the single-species pair anni- 
hilation process: the total density p(t) — J2iPi(t) ~ t -1 . 
In stark contrast, in d — 1 it exhibits species segregation, 
and is characterized by 



pit) ~ t 



-«(«) 



(1) 



with a q dependent exponent given, within the approach 
presented below, by a(q) = (q — l)/2q. We note that for 
q — > oo, two particles of the same species meet with zero 
probability; the distinction between the different species 
then becomes irrelevant, and this model is equivalent to 
the A + A — ► reaction f§], with known a(oo) = 1/2. 

In order to set the stage for our arguments, we briefly 
summarize the physics of the single- and two-species an- 
nihilation processes. For A + A — » the mean- field rate 
equation pa = —^p\ with the solution pA(t) ~ 1/Xt pro- 
vides a valid description only for d > 2. For d < 2 
nearby reactant pairs quickly annihilate, leaving only 
well-separated particles, which in turn slows the density 
decay down. These arziii-correlations mimic a repulsion 
between the particles; in a field theory representation of 
the associated master equation |g{ they lead to a down- 
ward renormalization of A. As the diffusion propaga- 
tor remains unchanged, the perturbation series is readily 
summed to all orders; one finds p A {t) - t- d / 2 for d < 2 
and p A {t) ~ t 1 Ini a t the critical dimension d c = 2 0. 

For the two-species pair annihilation A + B — > the 
rate equations read pa/b = —XabPaPb- With equal 
initial densities p A (0) = Pb(0) they are solved again 
by PA/B(t) ~ I/*! witn Pa(0) > Pb(0), say, one ob- 
tains pB(t) ~ exp(— Xab[pa{0) — Ps(0)]i) for the minor- 
ity species, while the majority approaches a saturation 
density pa(oo). In order to establish the effects of spa- 
tial fluctuations, it is crucial to notice that the density 
difference pa~Pb remains conserved under the reactions; 
for Da = Db it simply obeys the diffusion equation [[l0| . 
Initial density difference fluctuations therefore amplify 
in time relative to the total density. As a consequence, 
when pa(Q) = Pb{0), then for d < d c = 4 the system 
segregates into domains and the density decay is slowed 
down to Pa/b^) ~ t~ d / 4 jnj. The renormalization group 
provides a firm basis || for these arguments, at least for 
2 < d < 4: An effective field theory can be derived that 
corresponds to the rate equations plus diffusion terms, 
which establishes the segregation. 

For unequal diffusion constants, Da ^ Db, this pic- 
ture is not qualitatively altered; however, special initial 
conditions may change it. Consider, e.g., particles that 
initially alternate in one-dimensional space, ABAB..., 
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and that upon encounter react with probability one. The 
distinction between A and B is then meaningless and the 
system is in the A + A — > universality class For 
unbalanced initial conditions, pa(0) > ps(0), stretched 
]l2| rather than simple exponential relaxation ensues for 
rf < 2: lnp s (t) ~ -t d/2 , whereas \ap B (t) ~ -t/lnt at 
d c = 2 [0|. If the two species are already segregated 
initially, d c = 2 also is the sole critical dimension |ll| . 

This summary helps us classify the possible scenar- 
ios for the g-species MAM with arbitrary parameters 
Xij, Di, and Pi(0). Generically we expect that after some 
crossover time only the least reactive and initially most 
numerous species will have survived, resulting in an effec- 
tive two-species problem. After this reduction of q to the 
effective value g e ff = 2 the final asymptotic decay laws 
will be those of the two-species system with unequal ini- 
tial densities discussed above. However, on special sub- 
manifolds of parameter space, and in particular in the 
presence of symmetries, reduction to q e s = 2 may not 
be possible and novel behavior may appear. That not 
all symmetries lead to new behavior may be illustrated 
by the cyclic reaction scheme A + B — ► 0, B + C — > 0, 
C + D — > 0, and D + A — > 0, all with equal rates and ini- 
tial densities. Here, we may readily identify the species 
A = C and B = D, respectively, which takes us back 
to the A + B — > reaction with Pa{0) = Ab(0). In this 
Letter we address the most prominent case that requires 
special consideration, and will in fact lead to novel ef- 
fects, namely the fully symmetric MAM, in which all q 
species are equivalent (whence q e g = q). 

First, we notice that the renormalization of the annihi- 
lation vertices in this g-species model is independent of q 
and identical to that of the A + A — > reaction Q , with 
d c = 2 0]. Second, for q > 2 there exists no microscopic, 
local conservation law. As a consequence, following the 
arguments in Ref. ||, any memory of the initial state 
will eventually become lost. This eliminates the segre- 
gation mechanism at work in the q = 2 case. Next we 
invoke the mean-field rate equations and conclude that 
Pi(t)=p(t)/q~l/tford>2 0. 

For d < 2, however, one needs to extract the correct 
asymptotic scaling from the Callan-Symanzik renormal- 
ization group equation. This requires an explicit com- 
putation of the density (a function of its initial value) 
as a power series in the renormalized annihilation rate 
Xr. At the critical dimension d c = 2 the renormal- 
ized rate flows to zero logarithmically, A j? (t) ~ 1/lni, 
leaving merely the tree diagram contributions that cor- 
respond to the solution of the rate equation. Thus we 
predict p»(t) = p(t)/q ~ t _1 lni in d = 2. The diffi- 
culty for d < 2 is to demonstrate that for large values 
of the relevant operator p(0) the power series remains 
properly controlled |Q. For d < 2 this has proven elu- 
sive even for the two-species system Moreover, in one 
space dimension blockage effects of (hard-core) particles 
have been found to often play a crucial role in multi- 
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FIG. 1. Monte Carlo results for the total density decay 
vs. time in the pair annihilation reactions A + A — > and 
Ai + Aj -> with 1 < i < j < q (q = 2, 3, 4) for equal 
initial densities in two dimensions. The plots for q 7^ 2 depict 
p(t)/ln£; the solid lines indicate the functions t~~ x and t" 1 ^ 2 . 

species reaction-diffusion processes Jl5| . We shall see that 
the g-species MAM, too, develops entirely novel features 
when restricted to a linear chain: The q species segre- 
gate into well-defined domains, which remain stable be- 
cause of the mutual annihilation processes that prevent 
species mixing and the special one-dimensional topolog- 
ical constraints that do not allow a given species to in- 
teract with all others. As a consequence, we find that 
for all 2 < q < oo that the total density decays in d = 1 
according to the power law ([!]). 

We now present our numerical evidence in d = 1 and 
d = 2, and then proceed with the analysis of the one- 
dimensional model. In order to check the predicted uni- 
versal decay law in two dimensions, we performed Monte 
Carlo simulations on a 512 x 512 square lattice with hard- 
core particles. Starting from a full lattice with random 
distribution of q equally abundant species (q — 2,3,4), 
we let the particles perform unbiased random walks and 
annihilate with probability one upon encounter with a 
different species. One Monte Carlo time step was consid- 
ered complete after N(t) trials, with N(i) the number of 
remaining particles at that instant. The results, shown 
in Fig. |l|, clearly support p(t) ~ i -1 \nt for the g-species 
MAM with q = 3,4; this is similar to the decay law 
of the A + A — > reaction, and in contradistinction to 
p(t) ~ t^ 1 / 2 for q = 2. We have also checked the purely 
mean- field behavior for q = 3, 4 in a 50 3 cubic lattice p6[ . 

In d = 1, however, simulations of the MAM with equal 
initial densities yield decay laws that differ importantly 
from both the A + A -> and A + B -> cases [§. Fig- 
ure H shows our Monte Carlo results for q = 2, 3, 4, 5 on 
a chain of 10 sites with periodic boundary conditions. 
Evidently at long times p(t) ~ t~ at - q \ where a(q) in- 
creases with q from a(2) = 1/4 (the A + B — > value) to 
a(oo) = 1/2 (the A + A -> value). 

In order to study the one-dimensional MAM we sim- 
plify it such as to retain only its barest essentials. This 
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FIG. 2. Monte Carlo results for the total density decay 
p(t) vs. t in the pair annihilation reactions A + A — » and 
Ai + Aj -> with 1 < i < j < q (q = 2, 3,4, 5) and equal 
initial densities in one dimension. The solid lines indicate the 
functions t~ 1/2 and f" 1/4 . 

simplified version, to be referred to as SMAM, arises from 
the following considerations. The one-dimensional sys- 
tem may at any time be decomposed into a sequence of 
domains, each containing only a single particle species. 
Owing to the diffusive nature of the process the typical 
domain size increases as L(t) ~ {Dt) 1 / 2 . Let us assume 
the asymptotic decay law p{t) ~ p (p 2 Dt)~ a , where a re- 
mains to be determined. The average particle number in 
a domain then scales as n(t) — L(t)p(t) ~ [p^Dt)^ a+1 / 2 , 
and phase segregation occurs only if a < 1/2. Adjacent 
domains are separated by reaction zones, of which there 
are l/L(t) per unit of length. Therefore, as argued in 
Ref. jL7| for the two-species case, the total particle den- 
sity decreases as p(t) = — n(t)/L(t), with the typical 
number of annihilations per unit of time in a zone. The 
SMAM is now defined by the assertion that fluctuations 
in the annihilation rate n(t), whether in the course of 
time or between different reaction zones, are irrelevant 
and may be ignored; i.e., the particle content of each do- 
main, owing to the annihilations taking place at both of 
its ends, decreases at the uniform rate 2n(t). To com- 
plete the picture, we need to specify what happens when 
a domain loses all its particles: Then, with probability 
l/(q — 1), the left and right neighboring domains con- 
tain identical species, and consequently fuse into a sin- 
gle new domain; or, with the complementary probability 
(q— 2)/(q— 1), they contain different particle species and 
a new reaction zone appears. 

The SMAM may be cast into an explicit algorithm. 
We consider a one-dimensional lattice whose sites 1 < 
i < represent the domains of the original MAM. We 
randomly select initial values for the particle num- 
bers in each domain. This random initial state evolves in 
time via deterministic iterations. The (k + l)th iteration 
changes the total number of sites from N^ k ' to N^ k+1 ^ 
and converts the integer set {n- }!jt-i to 
by successive application of the following four iteration 



steps: (i) All nf^ are reduced by 1. (ii) All sites that as a 
result have become empty, are eliminated and the other 
sites are reconnected without reordering, (iii) Any two 
sites that as a result have become neighbors, are, with 
probability l/(q—l), fused into a single site whose num- 
ber variable is the sum of the number variables of the 
fusing sites, (iv) The remaining sites are relabeled with 
an index 1 < j < N^- k+1 \ The fcth iteration yields the 
total number of domains and the average particle 

number rffl per domain. The particle density and the 
physical time follow from p^/p^ = N^n^/N^n^, 
and i( fe )/t(°) = [N^/N^] 2 , as N(t) ~ L(t)- 1 ~ t' 1 / 2 . 
A key feature of the SMAM is that at every iteration 

(k) 

step k the numbers n\ are uncorrelated, for they de- 
scend from disjoint sets of 'ancestor' variables. Therefore 
this model obeys an exact closed set of equations, which 
we will now derive and analyze. Let us for the moment 
suppress the iteration superscript (k), and denote, pre- 
ceding the /cth iteration, the total number of domains by 
N, the total number of domains containing n particles 
by M„, and their relative abundance by f n = M n /N. 
Primed symbols will indicate the corresponding quan- 
tities after the fcth iteration. In one iteration the to- 
tal number N of domains diminishes by M% due to step 
(i). Step (iii) results in the additional disappearance of 
domains; calculating their exact number requires taking 
into account all instances where two or more vanishing 
domains form a sequence of nearest or next-nearest neigh- 
bors. After some combinatorics one finds |16| that N' 
and N are related by N' = (1 - /i) [1 - f t /{q - 1)]N. 
By means of more elaborate combinatorial analysis one 
may express the final number M' n of n particle domains 
in terms of the initial M m . The intensive variables /„ 
(n — 1, 2, . . .) then obey the recursion relation 



/; = [1 + (q - 2)A] (/„ +1 + hllr, 
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with the abbreviations /i = fi/[(q— 1)(1 — /l)] and lZ n = 

Eoo ps— 2 v->oo r p r 

s=2 Jl L-m%\,...,m a = \ JlTii+l ■ ■ ■ Jm s + 1 °n,m 1 + . . .+m s ■ 

The term of index s represents the creation of a domain 
of n particles by simultaneous fusion of s domains. The 
fusions with s > 3 are clearly very model specific and 
one would expect the essential physics to be embodied 
already in the lowest-order nonlinearity. Indeed, by trun- 
cating Eq. (|^) after the s — 2 term one obtains an elegant 
Boltzmann-like equation; the mathematical analysis be- 
low is easier, however, if all terms are retained. 

To find a solution to Eq. (^) we substitute an expo- 
nential distribution /„ = e(l — e)™ . The recursion then 
yields a similar distribution, but with a new parameter 
e' related to e by e' = e[l — e/(q — 1)]. For this solution 
fi = e, which may be substituted in the relation linking 
N' to N. Since n = 1/e, the total particle density obeys 
p'/p = N'e/Ne' = 1 — e. After restoring the iteration 
indices we obtain the pair of recursion relations 
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to be solved with initial condition < e^ -* 
for a random initial distribution e^ -* = (q - 



n(0) 



The solution of Eqs. 



< 1 [e-9; 
l)/q] and 

J) and (Q) determines and 
t(fc) = ^0)^(0) e (0)/ p (*) e (k)]5. the desired function p(t) is 

then obtained by eliminating the index k. 

We have been able to carry this through explicitly only 
in an asymptotic expansion for large k. Whereas its lead- 
ing order is readily evaluated, more attention is required 
to deal with the subleading correction. By analyzing the 
recursion relation (0) one finds that 
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where c is a function of e^/(q — 1). Analyzing Eq. (ji|) 
with (||) inserted then yields 
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with A = linifc^oo fc^ 1 IlLoI 1 _ e ^]- Expressing as 
a function of k and inverting leads to k(t) ~ (Ci) 1 / 29 — 
pog(Ct)+2 g logc-(g-l)(g-2)]/2 g withC = ^(q-l) 2 / 
e (o) t (o)_ Finally, upon substitution in Eq. (||), 
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with a{q) = 1/2 — l/2q < 1/2, confirming species segre- 
gation self-consistently and establishing Eq. (|l|) for the 
leading density decay of the SMAM. For q — 2 we re- 
cover a(2) = 1/4, and the limit q — > oo gives correctly 
a(oo) = 1/2. Notice that the term with log(Ct) and 
the dependence on c have canceled out in Eq. ([?]). The 
next-to-leading behavior is identical with the power law 
decay for the A + A — > reaction in d = 1. Its relative 
amplitude increases with g; thus it becomes numerically 
difficult to distinguish it from the leading term. We can- 
not establish that the correction term in Eq. (|7|) has the 
same relevance for the original MAM as we believe the 
leading-order term does; in the accessible time window of 
the MAM simulations our data are best described by an 
effective exponent a e g, which reflects a sizeable next-to- 
leading correction |l6j. Current large-scale simulations 
by Ben-Avraham and Zhong indeed confirm unambigu- 
ously both our leading density decay law (|J) as well as 
the power i -1 / 2 for the first correction term fll8| . 

In this one-dimensional system the reaction zone width 
£{t) is just equal to the typical interparticle distance be- 
tween representatives of different species. The reaction 
rate is just the inverse of the time needed to diffuse 
over this length [[l?]], hence n(t) ~ D/l{t) 2 . Combining 
this with p = —k/L and the known time dependences of 
L(t) and p(t), we find £(t) 2 - p Q 2 (pg£)i) Q! ^)+ 1 / 2 , whence 



with X(q) = (2q — l)/4g. For q = 2 this reproduces the 
known result A(2) = 3/8 @- The value A(oo) = 1/2 in- 
dicates that for infinitely many species the reaction zones 
grow as fast as the typical domain size, and hence there 
can be no segregation. How to aptly take into account 
the special one-dimensional topological restrictions in a 
field-theoretic description remains an open issue. 
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